AGN obscuration from winds: from dusty infrared-driven to 

warm and X-ray photoionized 



A. Dorodnitsyn 1,2 , T. Kallman 1 



ABSTRACT 



o 

u 



(N 

o 

> 
O 

We present calculations of AGN winds at ~parsec scales, along with the associated 
obscuration. We take into account the pressure of infrared radiation on dust grains and the 
interaction of X-rays from a central black hole with hot and cold plasma. Infrared radiation 
(IR) is incorporated in radiation-hydrodynamic simulations adopting the flux-limited diffu- 
sion approximation. We find that in the range of X-ray luminosities L=0.05 - 0.6L e dd, the 
Compton-thick part of the flow (aka torus) has an opening angle of approximately 72° — 75° 
regardless of the luminosity. At L > 0.1 L c && the outflowing dusty wind provides the obscu- 
ration with IR pressure playing a major role. The global flow consists of two phases: the cold 
flow at inclinations 9 > 70° and a hot, ionized wind of lower density at lower inclinations. 
The dynamical pressure of the hot wind is important in shaping the denser IR supported 
flow. At luminosities <0.1L edd episodes of outflow are followed by extended periods when 
the wind switches to slow accretion. 
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1. Introduction 

Active galactic nuclei (AGNs) are among the most fascinating objects in the universe: 
various energetic processes interplay in a volume a few cubic parsecs around a supermassive 
black hole (BH); radiation is likely decisive in connecting the accretion and feedback modes 
of AGN. Dusty, radiation driven flows are very effective in removing large masses of accreting 
gas, returning it back the galaxy, and naturally limiting the growth rate of a supermassive BH 
by the feedback with the host galaxy. Radiation also provides a non-linear causal connection 
between distant parts of an AGN itself, coupling together vastly different physical scales. 
Despite the pivotal role of radiation in the first principle modeling of AGN and the significant 
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progress in theoretical modeling over the last 10 years, proper incorporation of radiation into 
numerical models still presents a significant computational challenge. 

Well established phenomenology places AGNs into two distinct classes which roughly 
differ by the presence (type I) or absence (type II) of the broad emission lines in the optical 
and UV spectra. According to AGN unification paradigm ( Rowan- Robinson|1977 Antonucci 



& Miller 1985 Urry & Padovani 1995), this dichotomy is an artifact of a pure geometrical 



origin. In order to be observed as type II object, an AGN must be viewed at high inclination, 
in extreme cases close to the equatorial plane, while in the case of a type I AGN the observer's 
line of sight is at lower inclination, i.e. closer to the symmetry axis. The presence of 
geometrically thick toroidal obscuration is pivotal to the unification scheme of AGNs. 

The discovery of broad permitted emission lines in the polarized spectrum of the nearby 



Seyfert II galaxy NGC 1068 Antonucci (1984), and Antonucci & Miller (1985), provides 



evidence that a type I nucleus is harbored inside the obscuring ring of matter. The observed 
polarization is likely due to scattering of light emitted by the type-I core and then being 
scattered within an outflow emanating from the same nucleus. 

Observations of AGNs and quasars reveal that dust in large quantities can be present 
within a few parsecs from the BH, and a significant part of X-rays and UV are reprocessed 
into infrared (IR). For example, in luminous quasars almost half of the bolometric luminos- 
ity can be reprocessed into IR. This results in a 1-100/iin "hump" in the quasar spectral 



energy distribution (SED) Richards et al. (2006). In nearby AGNs multi-temperature dust 
distributions are observed directly by the Very Large Telescope Interferometer (VLTI). Such 



mid-infrared observations of the prototypical Seyfert II galaxy NGC 1068 (Jaffe et al. 2004) 



and the Circinus galaxy (Tristram et al. 2007) depicts a picture of a dusty, clumpy plasma 



exposed to external illumination. 

Various models have been proposed to explain the real or apparent thickness of the 
torus: 



1) a "warped disk model" , these are global warps in a locally geometrically thin accretion 



disk 



i.e.. 



Phinney 1989; Sanders et al. 1989) 



2) i) models of a quasi-static, rotating torus involving "clouds" (i.e., Krolik & Begel 



man 



1988 Beckert & Duschl 2004; Nenkova et al. 2008). Vertical thickness of the torus is 



supported via large velocity dispersion of clouds, and small-scale magnetic field is required 
to maintain elasticity of clouds during collisions; ii) models which rely on infrared radiation 



pressure on dust to keep rotating torus vertically thick (Krolik 2007). Most of IR radiation 



in this models comes from reprocessing of external X-ray, and UV radiation. 
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3) a "wind torus" A collar of outflowing matter may provide obscuration, i) MHD 
wind: possibility is that a global magnetic field can be involved in the form of an MHD wind 
( Konigl fc Kartje||1994 Elitzur fc Shlosman||2006 ) or directly supporting a quasi-static torus 



(Lovelace et al. 1998), or in a combination of radiative and magnetic driving (Emmering 



et al. 1992 Everett 2005; Keating et al. 2012); note that sufficiently strong (of the order 
of the equipartition one) poloidal magnetic field can also help to solve the problem of how 
the angular momentum is removed from accreting matter, however the origin of such highly 
ordered field remains unclear ii) an outflow (or failed wind) driven by IR radiation pressure 



on dust (Dorodnitsyn et al. 2011 2012). All these models have their benefits and drawbacks. 



Models investigating the properties of a self-gravitating, non-radiative, dissipationless torus 



through N-body simulations have been made by Bannikova et al. (2012). 



The accreting dusty gas at ~parsec distances from the BH is susceptible to self-gravitating 
instability; since it is exposed to intense radiation strong gravito-thermal instabilities can 
develop. 

In the outflowing wind, the non-linear phase of thermal instability is associated with 
"clumps" (i.e., Elitzur 2008). Presumably such clouds are confined by small scale magnetic 



field. Whether they are secondary to the flow and behave more like "waves", i.e. temporary 
structures evolving on a dynamical or radiation times scales, or real long-lasting "clumps" 
remains unresolved. All models involving clouds need to explain how to avoid too much 
dissipation through the cloud-cloud collisions, which would otherwise raise the temperature 
to a fraction of the virial temperature r vir g . In the present studies we do not consider the 
effect of magnetic field (both ordered and chaotic). 

It has long been recognized that the geometrical thickness of the obscuration presents 
a significant problem. In order to explain why the torus is geometrically thick, models have 
been constructed which rely on gas pressure, magnetic forces, and on radiation pressure 
either alone or in combination with other mechanisms. If one relies on gas pressure to 
support vertical thickness of the torus then virial arguments suggest that the temperature 
of such gas should be of the order of T vir g = 2.6 x 10 6 M 7 /r pc K, where M 7 is the black hole 
(BH) mass in 1O 7 M , and r pc is the distance in parsecs. Such high temperatures are not 
compatible with the existence of dust. 

As mentioned before, the presence of a small scale magnetic field can provide necessary 
elasticity to cloud-cloud collisions. However, this argument is not without a flaw. It is 
natural to assume that turbulence develops. If so, this leads to the magnetic field being of 
the order of the equipartition field. A result is that the dissipation of magneto-sonic waves 
will raise the electron temperature to be of the order of T v i r;g . 
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The connection between a toroidal obscuration, UV absorption and warm absorbers can 
be self-consistently assessed by global multidimensional simulations. The prequisite for such 
an endeavor is to include various physical effects and couple them to radiation. Additionally 
why such calculations are currently beyond the reach is that high numerical resolution is 
needed to resolve multiphase behavior of the outflowing gas. These are the reasons why 
different works including present studies assess smaller aspects of the bigger picture. 

The connection between a Broad Line Region and dusty outflow was considered by 



Czerny & Hryniewicz (2011). They suggested that the gas of the dusty accretion disk is 



weekly coupled to the BH due to radiation pressure on dust resulting in a dusty wind which 



gives the low-ionization part of the broad line region. Bottorff et al. (2000) considered the 



relation between the observed properties of UV and X-ray absorbers with those predicted 
from a magnetohydrodynamical self-similar wind model. They conclude that in case of the 
well-studied Seyfert 1 galaxy NGC 5548 the intrinsically clumped absorber is biased toward 
smaller radii in case of warm absorber and are likely located at higher distances in case of 
UV absorption features. 

Other models are based on the possibility that the pressure of the infrared radiation on 
dust may itself be enough to support the vertical thickness of the torus. The idea of a static 



infrared support was explored in Krolik (2007) where an approximate analytical model for 



the torus was constructed showing that there is enough radiation force to support a static 



vertically thick torus. However, it was argued in Dorodnitsyn et al. (2011 ) (Paper I), that it is 



very unlikely that the IR supported torus would be static. It was shown that the equilibrium 
between rotational, gravitational and radiation forces cannot be maintained if the temper- 



ature in the torus exceeds that of T vir r ~ 



-1U/4 



987(n/10 7 M 7 r^) 1 ^K, 
where n is the number density, necessarily resulting in the torus being an outflow. In the 



312 (n/10 5 M 7 r 



pc 



following paper (Dorodnitsyn et al. 2012 Paper II) we supported this idea with radiation- 



hydrodynamics simulations of the outflowing accretion disk wind, driven by the radiation 
pressure on dust. The distribution of IR radiation was calculated using a flux-limited diffu- 
sion approximation. 

An important simplification was made in Papers I and II: Only the infrared- driven 
part of the wind was considered. The argument was that this portion of the wind is far 
enough from the BH so all the conversion of X-rays to IR already happened outside the 
computational domain. That allowed us to consider only one group of radiation (the IR). 
The radiation temperature at the boundary was calculated to match the input energetics 
of soft X-rays. In other words we did not treat explicitly the reprocessing of soft X-rays 
into IR. Correspondingly, we were not able to follow the interaction of the X-ray heated hot 
component of the flow with the cold, dusty IR supported component. 
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In this paper we calculate the wind structure and explicitly treat the interaction between 
X-rays and gas with subsequent conversion of X-rays to IR. This allows us to follow in 
dynamics the transition from hot, photoionized to cold, dusty and IR supported flow. To 
make calculations numerically tractable we made several major simplifications in treatment 
of the transfer of X-rays as well as their interaction with the gas. For example, transfer 
of X-rays is reduced to simple attenuation with no multiple scatterings taken into account, 
and we approximately calculate the heating of dust grains by X-rays. The rest of our 
computational radiation-hydrodynamics framework remains unchanged from Paper II: the 
transfer of the infrared radiation and its interaction with matter is treated in a flux-limited 
diffusion approximation. The calculations are time-dependent and 2.5D, i.e. 3 dimensions 
in cylindrical coordinates and assuming axial symmetry. 

The plan of this paper is as follows: After this introduction the flow of this paper goes 
through Section 2 where we review basic assumptions of our model, including the description 
of the radiating fluid of the AGN dusty wind by means of radiation hydrodynamics equations. 
Interaction of the cold and hot components of the wind with X-rays are described in Section 
3 including corresponding approximations. In Section 4 we describe the numerical setup, 
initial and boundary conditions. The results are summarized in Section 5 while the final 
discussion is postponed to Conclusions, where we highlight prospects and limitations of our 
approach to the problem of AGN winds. 



2. Assumptions and physical model 



The flow in which the continuum radiation plays an important role is described by the 
system of radiation hydrodynamics equations. To first order in v/c, these can be formulated 
in the following form (Mihalas & Mihalas 1984): 



D t p + pV-v = 0, (1) 
Av = — Vp + g rad -W, (2) 

P D t(j\ = -pV -v -4n XP B + c XE E + n 2 (T -A), (3) 

pD t (j\ = - V • F - Vv : P + 4n XP B - c X eE, (4) 

where quantities related to matter: p, p, and e are the material mass density, gas pressure 
and gas energy density, v is the velocity; quantities related to radiation are the frequency- 
integrated moments (for definitions see Appendix A): E is the radiation energy density, F 
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is the radiation flux, and P is the radiation pressure tensor; Xp> Xe are the Planck mean 
and energy mean absorption opacities (in cm -1 ), c is the speed of light, B = aT A /ir is the 
Planck function, and a = ac/4 is the Stefan-Boltzmann constant, T is the gas temperature; 
other notation include the convective derivative, D t = + v • V, and Vv : P denotes the 
contraction (djV^P^ . Notice that all the dependent variables in ([lJ-Q are evaluated in the 
co-moving frame. Notice that E = aT A , where T r is the radiation temperature. Interaction 
between X-rays and the gas is incorporated in the last term of the gas energy equation (J3]), 
where F and A are respectively heating and cooling rates. The equation of state for the gas 
is assumed to be polytropic: p = Kpj, where 7 = 1 + 1/n, and n is the polytrope index, 
and p = (7 — l)e. The gas is one-component, one-temperature T = pp/plZ, where /1 is the 
mean molecular weight per particle, 1Z = 8.31 x 10 7 ergK -1 g -1 is the universal gas constant 
and plasma with 7 = 5/3 is assumed to constitute the flow. Three components of the flow 
velocity Vr, v z , and v§ are calculated in cylindrical coordinates z, R, assuming azimuthal 
symmetry, = 0. 

Frequency- independent moments E, F, which appear in the above set of RHD equations 
are obtained by calculating angular moments from the frequency-integrated specific intensity, 
7(r, Q, v, t) and are given in Appendix A. Forces which are explicitly taken into account 
include radiation, g ra d pressure and gas pressure, p _1 Vp. 

These equations should be supplemented by the radiation force, g ra d, which is calculated 
from the following relation 



grad = , 5 

C p 

where Xf = Xa + Xt is the total flux mean opacity consisting of absorption opacity, Xa and 
the Thomson scattering opacity, xt- As in Papers I, and II we did not differentiate between 
Xf, Xp an d Xe- Radiation pressure on dust is the most important from the dynamical point 
of view (Paper I). Notice that in the infrared, the Rosseland mean opacity, K d = Xd/p of 
dust with the temperature 10 2 — 10 3 K, is approximately 10 — 30 times larger than that 



of the electron Thomson opacities (Semenov et al. 2003). In this paper when calculating 



radiation pressure force in ^ we take into account only IR pressure on dust. Including 
UV force (in lines and on dust) is left for future studies. We adopt a simple procedure for 
the calculation of g ra d: if the temperature of the dust, T d > T sub , where T sub is the dust 
sublimation temperature (see further in the text) the force is assumed to be zero. 

Our description takes into account two types of radiation: infrared and X-ray. 

The transfer of X-rays is treated in a single stream approximation. This implies that 
X-rays are simply traced from the point source (i.e. corona) adopting e~ Tx attenuation, 
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where Tx is an optical depth with respect to X-ray absorption ( Dorodnitsyn et al. 2008). 
The distribution of the infrared radiation is calculated in a Flux Limited Diffusion approxi- 
mation. Detailed discussion about the validity of equations ([lJ-Q is given in Paper II. We 
briefly remind that the closure relation between F and E is obtained adopting the diffusion 
approximation: 

F = -DVE, (6) 
where if optical depth is large r>l, the diffusion coefficient reads: 



D = c\, (7) 

where A = is the photon mean free path, and k — x/P- From ([5]), ([6]), and ([7]) one 

can see that in the diffusion regime the radiation force does not explicitly depend on the 
opacity. 

The diffusion approximation tacitly assumes that optical depth r > 1. Most of the 
torus where IR pressure is important indeed has Td > 1 where Td is the optical depth of the 
gas-dust mixture in the infrared. 

If r < 1 the diffusion approximation should be modified so a correct limiting behavior 
at t < 1 is regained. One can see that without such modification when r C 1 the mean free 
path, A — > oo, and D — > oo, and |F| — > oo which is in contradiction with a free-streaming 
limit where it should be |F| — > cE ([7]). That is, when optical depth becomes small, or when 
p — > 0, the standard diffusion approximation is no longer applicable. 

In order to describe correctly regions of small r, the standard approach is to adopt 



the flux-limited diffusion approximation (Alme & Wilson 1974 Minerbo 1978 Levermore 



& Pomraning 1981). In this approximation A is replaced by A* = AAp, where Ap is the 



flux limiter. The flux limiter we adopted in Paper I, II and in the current work is that of 



Levermore & Pomraning (1981): 



A F 



2 + Rlp 



6 + 3i? L p + R 2 



LP 



where i?LP = A \VE\/E. One can see, that if r — > 0, then i?LP — > oo, and \F\ 
t > 1 R LP ->• and A F ->• 1/3. 



cE\ and if 
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3. Interaction of radiation and the wind 

3.1. Photoionization equilibrium 

The premise of the current work is that the deposition of energy at parsec scale is 
dominated by external sources, i.e. by UV and X-rays from the central black hole and inner 
accretion flow. Due to the high opacity of the dust-gas mixture, a significant part of the 
UV and soft X-rays is absorbed and reprocessed into infrared in a thin layer of thickness, 
5l/R lpc ~ 1.3 x KT 3 ^ 1 (Paper I). 

The dynamical time within the flow is usually much larger than the characteristic time 
of the photoionization and recombination. Thus, the ionization balance is determined by the 
condition of photo-ionization equilibrium. The condition of a photo-ionization equilibrium 
in the wind implies that the state of the gas can be be reasonably parameterized in terms 
of one "ionization parameter", £ -the ratio of radiation energy density to baryon density 



(Tarter et al. 1969). In a popular form it reads 



£ = 4 7T Fjn ~ 4 • 10 2 • / x T M 6 /(iV 23 r pc ) , (9) 

where F x is the local X-ray flux, L x is the X-ray luminosity of the nucleus, and n is the 
gas number density; N 2 s is the column density in 10 23 cm -2 , f x is a fraction of the total 
accretion luminosity Lbh available in X-rays, T is a fraction of the total Eddington luminosity 
-^edd = 1.25- 10 44 Mq, where Mq is the mass of a BH in 1O 6 M . We adopt a simple attenuation 
model in which the local deposition of energy is proportional to F x = L x e~ Tx /(47rr 2 ), where 
the optical depth at soft X-rays is integrated over a straight line connecting the point source 

corona with a given parcel of gas r x = / K x pdl, and k x is the X-ray opacity. 



3.2. Hot and cold gas 



To calculate photoionization balance, we take into account Compton and photo-ionization 
heating and Compton, radiative recombination, bremsstrahlung and line cooling. The rates 



of these processes are calculated making use of the XSTAR photo-ionization code (Kallman 



fc Bautista||2001 ) for the incident spectrum which is a power law with energy index a. As the 
direct incorporation of the XSTAR subroutines into the hydrodynamics code is not feasible, 
the key ingredient of our method is to make use of approximate analytic formulas for the 
heating and cooling rates. These formulae approximate extensive tables of heating-cooling 



rates obtained from XSTAR; these approximations are similar to those of Blondin (1994) 



although modified by Dorodnitsyn et al. (2008) to incorporate better atomic data and power 
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law ionizing continuum with a—1: 

^(ergcm^^ 1 ) = 8.9 • KT 36 £ (T x - 4T), (10) 
for the Compton heating - cooling; 

r x (erg cm 3 s" 1 ) = 1.5 • lO" 21 ^ 1/4 T' 1 ' 2 ^ - T)T~\ (11) 

for the photo-ionization heating-recombination cooling, and for the bremsstrahlung and line 
cooling: 



A (erg cm 3 s 



3.3 ■ 1(T 27 T 1/2 

(4.6 • 1(T 17 exp(-1.3 • i 5 /T)£ ( ~ a8 -°' 98Q) T-l/2 + 1(T 24 ) 5. (12) 



Formulae (10)-(12) were found to be in reasonable (25%) agreement with numerical results. 
A variety of physical and chemical processes influence the transformation of X-rays 



into IR if dense dusty, molecular gas is illuminated by external X-rays (Maloney et al. 



1996 Krolik & Lepp 1989). Rather than attempt to incorporate all these processes in our 



radiation hydrodynamics framework, in the current paper we adopt several approximations 
which reproduce the qualitative behavior of low temperature gas exposed to external X- 
ray radiation. At column densities, A" smaller than 10 24 cm -2 photons with energies below 
10 keV interact with the gas through photoelectric absorption while Compton scattering 
dominates at higher energies. The total ionization per particle at a given point in the 
wind is proportional to £ e ff, the ionization parameter which takes into account the frequency 



dependent attenuation of X-ray flux. It is shown by Maloney et al. (1996) that heating 



and cooling rates of such X-ray illuminated molecular gas depend on the modified ionization 
parameter: £ e ff = ^/A^ 9 = 1.26 x lO^ 1 F x / (ngN®^) , and several regimes are approximately 
distinguished: 1) The regime of highly ionized gas £ e ff ^> 1. In this regime for the heating 
and cooling rates we make use of approximations (10) - (12); 2) If 10~ 3 <C £ e ff < 1 the gas 



is primarily atomic, and we also adopt heating-cooling rates (10) - ( |12| ); 3) If the ionization 
parameter, is small (£ e g- <C £ m = 10~ 3 ) the gas is largely molecular. This is where the most 
complexities due to chemistry and radiative interaction with the dust occur. The proper 
incorporation of this regime calls for consideration of all the radiative, chemical and dust 
network of processes. This is challenging even without coupling to hydrodynamics. 

In the regime when £ e ff < £ m , we assume that the molecular gas is in radiative equi- 
librium with dust, i.e. T g = where equilibrium dust temperature is found from the 
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attenuated X-ray flux, F x . Thus we approximately take into account that dust can directly 
reprocess X-rays to IR contributing to the IR energy density, E. We approximately take this 
contribution into account, using equation (see, Appendix B): 



dT r 
dt 



n d c a d (T d - T r ) 



(13) 



Dust temperature is found from the approximate relation, Td = (4Tx/(ac)) 1,/4 (notice that 
E = aT r 4 ), where Fx is the local attenuated X-ray flux. If the dust temperature Td > T su b, 
where T sub is the dust sublimation temperature then it is assumed that dust is destroyed 
and no conversion to IR occurs. 

In our approach, the main difficulty in calculating the energy density of IR radiation 
from X-rays resulting from a single fluid approximation. By having only one species (gas) 
and not evolving a separate dust component, we can mimic the interaction between dust 



and X-rays, and calculate Tj by means of equation (13), while additionally assuming the 



relation between Tj, and T in a dusty-molecular phase, in which, for simplicity we assume 
T = T r , with the later calculated from(13). Other dust characteristics in (13) are as follows: 
nd = p/d/^d is the dust grain number density, where fd is the dust-to-gas mass ratio, and ma 
is the dust grain mass, and a d = ixr\ is the dust cross-section, and r d is the dust grain radius. 
Throughout our calculations, we adopt the following fiducial values for the parameters of 
dust: typical dust grain radius, rd = 1 x 10~ 4 cm, the density of the dust grain pd = 1 g cm -3 , 
fd = 1 x 10~ 2 , and above the sublimation temperature T su b = 1500 K the dust is assumed to 



be destroyed. Equation (13) for Td(Tx) approximately describes how X-rays are converted 
into IR through reprocessing by dust. 

If £ c ff ^ £m and the T < T su b we assume that the exchange between T X (E) and T(e) is 
described by the balance between D t (E/p) and the last two terms in equation Q. The gas 
temperature T can contribute to IR and vice versa. For example, gas can contribute to the 
X-ray-IR conversion through cooling. Cooling can be due to radiative processes, i.e. the low 



energy limit of ( 11 )-( 12 ) and due to "adiabatic" cooling due to fluid motion. 



The update of E from (13) is performed adopting operator splitting in the source step, 



i.e., before the update of the equation Q is done. During the update over dt, the new 



E t (t + dt) is calculated from (13). 
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Methods 



To solve the system of equations ([T])-((4]) we expanded the radiation hydrodynamic frame- 
work described in our previous papers. The additions include the interaction of X-rays with 
dust grains and with gas. The coupling between the X-ray heating and cooling with hydro- 



dynamics has been implemented and tested in our earlier works (Dorodnitsyn et al. 2008). 



Here we have added physics which includes the direct interaction of X-rays with dust grains. 



The radiation-hydrodynamic equations are solved using the ZEUS code (Dorodnit 



syn et al. 2012) with extensions taken from Dorodnitsyn et al. (2008). The radiation- 



hydrodynamic part is designed and tested in (Dorodnitsyn et al. 2011, 2012) and conforms 



with methods and structure of the original ZEUS code (Stone & Norman 1992). 



In the following we briefly review the organizational structure of the code. Similarly 
to the way pure hydrodynamics is treated in ZEUS, the radiation part is also split into the 
source and transport steps. Artificial viscosity is adopted to resolve discontinuities in the 
flow. In the source step the difference equations for dE/dt and de/dt are solved, and in the 
transport step the advection of E and e is calculated. The update of E and e is calculated 
adopting fully implicit scheme (see Paper II); this includes the solution of the diffusion 
problem and the e -H- E update i.e. the update of gas and radiation energy densities due to 
radiation-gas interaction. 

Similarly, in the X-ray-matter interaction step, heating 7 X — > e, and cooling e — > 7 X 
due to Compton, radiative recombination, bremsstrahlung and line cooling is also done fully 
implicitly in a separate update. 

We make use of the R, z cylindrical grid, which extends from R[ n to R ont in radial and 
from z- in = to z out = R out in the vertical direction. The lower boundary at z- m = is 
intended to represent a thin accretion disk, which serves as an infinite reservoir of gas and 
dust. Wind can escape from all other boundaries, in principle, though the centrifugal barrier 
effectively prevents escape from the inner boundary. A fixed grid of size 200 x 200 is adopted 
for all calculations. Hydrodynamic boundary conditions (BC) are as follows: at the left 
(innermost), right, and upper boundaries we adopt outflowing boundary conditions. 

At the equator we fix the inflow speed v j at a small fraction of the speed of sound leaving 
the density to adjust according to the continuity equation. Hereafter, index j varies over R 
and index i over z. We allow p is j to vary and fix Vi S j, that is pi S j = Vi S+ ijPi S+ i t j/vi S , and Vi S = 
const, provided v is <C v is+ i :S , where v is+ i iB is the sound speed; everywhere in the simulations 
we assumed Vi S /vi S+ i yS = 0.01. Notice that this is different from the implementation of the 
equatorial BC in Paper II where it was assumed pi S j ~ RJ P and velocity was found from the 
continuity equation and only checked after the fact that Vi S j <C t>j S)S . Present implementation 
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of boundary conditions allows the disk to be flexible choosing mass loading of the wind. The 
initial distribution of density at the equator was taken to be Por' 1 , where p — m P n o is the 
characteristic density, m p is a proton mass. We find that numerical solution quickly forgets 
about initial equatorial distribution of p, and that solutions are not sensitive to the choice of 
p. Theoretically, since pi S j is allowed to change, the final solution should only weakly depend 
on initial equatorial distribution of p. Practically, however, this is true only for those parts 
of the disk where the wind is developed. We also find that the present implementation of BC 
produces a more stable solution near the equator, thus somewhat relaxing the restriction on 
the hydrodynamical time step. 

The radiative BC are free-streaming F ~ cE at all boundaries except at the equator 
where zero flux BC are implemented. 

5. Results 

5.1. Dynamics of the wind 

Deposition of radiation energy into the flow is controlled by the parameter r\ = Lboi/ -^Edd- 
We calculate models for rj = {0.05, 0.1, 0.3, 0.4, 0.5, 0.6}. The innermost boundary of the com- 
putational domain is located at radius, R = 0.5 pc. Other parameters are M B h = 1 x 1O 7 M , 
fx = 0.5, and the initial scale density Uq = 1 x 10 8 cm -3 . The evolution of the radiative flow 
was followed for 10-30 dynamical times, t ~ 1.5 x 10 11 r^M^^s. 

In Papers I and II the temperature of the IR radiation was prescribed at the boundary, in 
a physically motivated way (to match the energetics of X-rays, f x L ho \). The results obtained 
in Paper II are generally in a good agreement with the present modeling. However, explicitly 
including X-rays into calculations reveals several important features of the flow which were 
absent in previous studies. 

One of the results of the present study is that the pressure of the hot photo-ionized 
component is very important in shaping the IR flow. Excitation of such a hot, photoionized 
wind automatically mirrors the creation of the IR flow. The interface between the two 
components of the flow is seen in all simulations presented in this paper. In all cases it 
takes the form of oblique shocks or contact discontinuities. A supersonically moving hot gas 
obliquely hits the slowly moving IR flow creating a discontinuity (in most cases an oblique 
shock wave). 
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5.1.1. Model with L = O.GL e dd 

Figure [T] shows the development of the radiation-driven disk wind for T = 0.6. The 
outflow develops quickly after about one dynamical time, i ~ 4 x 10 4 yr. The mass-loss 
rate is (M) ~ 0.1 — O.2M yr _1 on average, and peaks at (M) ~ lM yr _1 in the period of 
time between 4 x 10 3 yr and 5 x 10 3 yr. The total computation spans from to 8 x 10 3 yr. 
These averaged peak numbers include both failed wind and wind which has enough energy 
to escape to infinity. 

The energetics of the wind is measured by computing the kinetic output of the wind 
at the outer boundary of the computational domain E: v^ n ~ 2L^ Q /M, where L^ n = 

/ pv 3 /2 dTi is the kinetic luminosity of the wind. For a radiation-driven wind we expect low 
Jt, 

values of Lkin/^boi; here we have Lkm/^boi — 5 x 10~ 5 — 2 x 10~ 4 . 

The averaged dynamics of the wind, is described by the the average bulk velocity of the 

flow (v) = / pvdV I / pdV, where V is the total or partial volume occupied by the flow. 
Jv Jv 

Here, the dense part of the flow has (v z ) ~ a few x 1km s _1 , and (vr) ~ a few x lOkms -1 . 
The maximum velocity reaches 1000 km s^ 1 in R, and over 600 km s" 1 in the z direction. 

Figure [2] shows a color plot of the distribution of the gas temperature. The gas spans 
a vast range of temperatures from being in a cold molecular state to photo-ionized, high 
temperature warm absorber flow. An apparent feature of our models is that the hot wind 
consists of large scale inhomogeneities. The cold flow doesn't show such large scale structure. 
Notice that we have a very simple test if the gas is in the cold, molecular-dusty phase, 
£eff < £m, and that if hot component is not shown the density structure is much more 
pronounced in the temperature plot. 

It is of interest to address the question of whether there can be the cold gas in a hot 
wind as well. However, the limited resolution of our studies does not allow us to provide a 
reliable answer to this question. We do not find the coexistence of hot and cold phases of gas 
in any appreciable quantities, but this does not exclude the possibility of such coexistence 
on length scales smaller than we can resolve. 



5.1.2. Model with L = 0.5L e dd-' the importance of the hot flow pressure 

Figure [3] which shows a color plot of the density at different times shares a lot of similarity 
with Figure [TJ One can distinguish approximately three regions in the flow: 1) a disk-like, 
geometrically thin and dense region, extending to R ~ 1.5 pc; 2) the IR-driven flow with 
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Fig. 1. — Color plot of the density, logp in (g cm -3 ) for the model with L = 0.6L ec id shown 
at different times given in years. Axes: horizontal: z: distance from equatorial plane in 
parsecs; R: distance from the BH in parsecs; 
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lower density, separated by contact discontinuity from 3) a hot, photoionized flow. As L 
decreases, so does the energetics of the wind. However, the vertical thickness of the dense 
IR supported flow is determined by the delicate balance between a) the pv 2 pressure of the 
hot component from above, and b) by the pressure of IR supported wind from below. 

Decreasing T from 0.6 to 0.5 causes the pressure of the hot flow to be less than in the 
solution shown in Figure [3] The balance between the pressure of the photo-evaporated wind 
and IR drive dusty flow is is slightly shifted towards IR-dominated flow. 

Figure [4] shows the 2D distributions of the velocity components v z and vr. The slowly 
outflowing wind is clearly seen; its well defined boundary extends to 0.5 pc. The flow consists 
of a "core", i.e. slow moving wind with radial velocities of a fewxlO < lOOkms -1 , and v z 
much smaller, of the order of the a fewxlkms -1 . Fast lower density component of the 
flow occupies most of the domain having velocities of the order of a fewxlOO < 700 km s -1 . 
The overall characteristics of the velocity field are similar to those of the model with L = 
0.6L edd . Not surprisingly, pumping less energy initially in X-rays generate less powerful 
wind: L kin /L hol ~ 1 x 10~ 5 - 7 x 10~ 5 . 

The mass-loss rate (M) ~ O.lA^yr^ 1 with the peak (M) ~ l.lM yr _1 in the period 
of time between 8 x 10 3 yr and 10 x 10 3 yr, and the computation spans from to 8 x 10 3 yr. 




Fig. 2. — Color plot of the gas temperature, logT in K for the model with L = 0.6L e dd, 
after 4.2 • 10 5 yr. Axes: horizontal: z: distance from equatorial plane in parsecs; R: distance 
from the BH in parsecs; 
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Fig. 3. — Color plot of the density, logp in (gem -3 ) for the model with L = 0.5L edd , shown 
at different times given in years. Axes: horizontal: z: distance from equatorial plane in 
parsecs; R: distance from the BH in parsecs; 
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Fig. 4. — The surface plot of the z— and R— velocity components. Model is shown with 
L = 0.5L e dd, t = 6.5 x 10 5 yr; Horizontal: R: distance from the BH in parsecs; z: distance 
from the equatorial plane in parsecs; 
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5.1.3. Model with L = 0.3L cdd 

Figure [5] demonstrates the evolution of the density for this model at different times. The 
disk-wind system goes through several stages: the domination of the hot wind; puffing up 
the IR supported disk; squeezing of the later towards the equator by the hot wind; building 
a high density disk-like outflow, etc. As in previous examples, the system doesn't come to a 
quasi-stationary state instead going through such episodes all over again. 

An interesting feature which is seen at t — 2.5 x 10 5 yr is in fact the residue from the 
contact discontinuity which existed at earlier times. The average velocity (v) ~ a few x 
172 km s -1 . This is about 60% of the escape velocity, U esc ; the numbers are given for t = 
20tdyn- The maximum velocity of the dusty flow is 242 km s -1 . 

The value of the mass- loss rate is remarkably similar to models with higher T: (M) ~ 
0.1 - O.2M yr" 1 reaching 1M yr" 1 . 

5.I.4. Models with L < 0.1L e dd- slow accretion with episodic outflow 

When the radiation input falls below an approximate threshold value no rigorous outflow 
is observed in our simulations. Instead we see that episodic outbursts of the hot evaporative 
flow are excited from the inner parts of the disk. Here we found this threshold luminosity to 
be L ~ 0.1L edd . The results for the model for L = 0.1L edd is shown in Figure [6j 

The density of this wind is too low to provide any considerable shielding from X-rays. 
Without much shielding it is difficult to form cold and dusty phase which can be accelerated 
by IR. Such an episodic outburst of the dense wind is followed by a gradual fall back of the 
dense component mostly in z direction towards equatorial disk. 

The velocity of the hot wind, v < U esc , and it skirts the denser flow and falls back towards 
equator at larger radii, R ~ 2.5pc. The density piles up there and the shielding from X-rays 
rapidly increases. Correspondingly, the ionization parameter drops below £ m and the cold 
IR-driven wind again has a chance to develop. Notice that the vertical component of the 
radiation force is roughly g e ff, 2 ~ (dE/dr) cos 6, where r is the spherical radius, and 9 is the 
inclination angle, measured from z axis. Since the vertical extent of the IR wind is small, 
g e s,z is small, and the dynamic pressure of the IR wind cannot overwhelm the pressure of 
the hot component at higher z. The IR wind stalls and starts to slowly fall back (accrete) 
towards the equator. Density outside this dense IR-supported "pancake" drops and the gas 
gets overionized again. The cold disk itself shrinks from both inside and outside. 

We also calculated a model with L = 0.05L edd , and found that the general behavior of 
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this model is similar to previous with L = 0.1L e( jd- The wind is episodic and intermittent 
with the "IR-supported pancake" phase. 

It is difficult to access the dynamical properties of models with L < 0.1L e dd, because 
our equatorial BC excludes the influx of matter into the accretion disk. Strictly speaking by 
assuming thin accretion disk as a source of matter for the wind together with outflowing BC 
conditions we assumed that the wind, even if very weak, is always present. We will further 
elaborate on this limitation in the following section. 




Fig. 5. — Color plot of the density, logp in (gem -3 ) for the model with L = 0.3L e dd shown 
at different times given in years. Axes: horizontal: z: distance from equatorial plane in 
parsecs; R: distance from the BH in parsecs; 
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Fig. 6. — Color plot of the density, logp in (g cm -3 ) for the model with L = 0.1L e dd shown 
at different times given in years. Axes: horizontal: z: distance from equatorial plane in 
parsecs; R: distance from the BH in parsecs; 
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6. Obscuring properties 

We calculate Thomson optical depth, tq = j, K e pdl, where I is measured along the line 
of sight at the angle, 9 from the vertical axis, and averaging is made over all calculated 
models for particular V. Processing different models, we are interested in an angle 8 p ^ where 
r ph = r (#ph) = 1? i-e. when the wind becomes opaque. 

In the following we summarize results, listing models in the order from low to high 
luminosities. When V = 0.05 the optical depth rises sharply from 0.1 to 0.5 at (6) ~ 75°, 
(hereafter # e dge) and the photosphere is at (6* p h) — 78°. 

Increasing the luminosity to T = 0.1 causes the dense part of the wind to become 
thicker in vertical extent, optical depth rises at the angle 8 e d ge — 73°, and the photosphere 
is at (0 ph ) ~ 75°. At T = 0.3 optical depth rises sharply from 0.1 at (9) ~ 73°, (0 ph ) ~ 75°. 
and ^edge — 73°. 

For the model with T = 0.5 the photosphere is at (0 p h) — 75° and the torus edge is at 

^cdgc — 73 . 

Figure [7] shows a color plot of a radial Thomson optical depth measured from a BH to 
a given point of the computational domain. The model shown has L = 0.6L e dd- This is the 
most luminous model we have calculated. It has the photosphere located at (9 p b) — 72° and 
the torus edge is at ^edge — 71°. 



Mass-loss rate 



The opacity of dusty plasma in the infrared domain is 10-20 times that of electron 
scattering (Semenov et al. 2003). Such huge opacity makes possible for the AGN to have 



a strong wind, radiating at a fraction of Eddington luminosity (defined with respect to 
Thomson scattering (see Paper I for an estimate of an outflow onset) 

Mass-loss rate from the IR driven portion of the wind can be estimated combining the 
results from a theory of stellar winds and some of the approximations adopted in Paper I. 
From the stellar wind theory the mass-loss rate from a spherically symmetric wind can be 



found from the following approximate relation (Lamers & Cassinelli 1999): 



M/(4tt) ~ Lm/c (r IR - l)/(4vrt; 00 r IR )) r w 



(14) 



where Tir = Lx/L c ad,m, v°° is the wind terminal velocity, and r w is the wind optical depth 
in IR. 



Fig. 7. — Color plot of the radial Thomson optical depth, log r from a BH towards a given 
point. Axes: horizontal: z: distance from equatorial plane in parsecs; R: distance from the 
BH in parsecs; 



For simplicity, as in Paper I in this section we assume that all incident UV and X-ray 
radiation is reprocessed within a narrow conversion layer which effective temperature, T e g can 
be found from the following approximate relation: aTrjxF e dd = crT^. where T = L/L ec jd, 
r]x = 0.5 is the fraction of X-ray radiation from the total radiation , and the fraction 
a ~ 0.5 of the incident flux is re-emitted outwardly. Adopting the above, one obtains 
Tir = ttT^x^d/^e which is approximately 1.25 for M B h = 1O 7 M , L = 0.5L ddd and K d = 
10n e . The value of v°° is estimated while neglecting P g in favor of the radiation pressure: 
v°° = ^/2GM/r (riR — 1) ~ 207 km s -1 adopting the above set of parameters plus the wind 
launching radius Ro = 0.5 pc, and assuming that the wind occupies a wedge-like (in a quarter 
of the domain) space of an opening angle of 75° we arrive to the estimate: M ~ 3 t w M yr _1 . 
Monte-Carlo simulations show that the factor r w can reach 1-5 for an AGN radiating at a 
fraction of L ed d (Roth et al. 2012), meaning the potential for the wind tori to reach mass 
loss rate of lOM yr _1 ; We do not see such mass-loss rates in our simulations, which may 
indicate that instead of being blown away as it would certainly do in a spherically-symmetric 
case, the gas is effectively ablated towards the equatorial plane by the pressure of the hot 
gas component. As the latter is itself the result of interaction with radiation, one can say 
that radiation effectively clears its way through the opening of the torus. 
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8. Discussion 

The following average properties of our wind solutions can be summarized: 

. The total mass-loss rate in the IR driven and warm absorber flows is between (M) ~ 
O.lM yr -1 and (M) ~ l.5M Q yi~ l depending on the time of the evolution of the 
wind. This includes also the failed wind. The amount of gas which has enough energy 
to escape to infinity is in a good accord with the energetics derived from accretion. 
Notice that approximately 1.6 x 10 _1 — 1.3 x 1CT 2 M Q yr _1 are required to be accreted 
through a geometrically thin accretion disk to support our models. Also there is no 
problem of depletion of the torus since we have an infinite reservoir of matter in a razor 
thin equatorial disk. 



2. The IR driven- wind has the velocity in the range of of a fewx 10 
fast hot wind has velocity in the range of 400 — 800kms _1 . 



200 km s l , and the 



The transmission properties are similar between different models. The characteristic 
angle at which the wind becomes opaque to Thomson scattering is about 72° — 75°, and 
is determined by the balance between the hot, photo-ionized wind and the cold wind 
supported by the IR pressure on dust. It is interesting that the geometry of the torus 



is close to one obtained in Dorodnitsyn et al. (2008) at late times when the pressure 



of the hot, evaporative component becomes important. 
4. Both hot and cold components of the wind are non-stationary. 



Models such 


as ' 


Konigl & Kartje 


( 


1994 


); 


(1992 


); 


Everett ( 


200, c 


>); 


Keating et al. 


(2012) 



Elitzur & Shlosman 


(2006) 


; Bottorff et al. 


(2000) are 



premise of 1) the existence of a strong global magnetic field, and 2) uses the prescription 
of self-similarity. To a certain extent this limits the value of the direct comparison of our 
results with theirs. Indirect comparison would include: performing the detailed analysis of 
the spectral properties of the IR-driven "wind torus"; performing our simulations for the 



parameter range attributed to a well studied source (such as NGC5548, in case of Bottorff 



et al. (2000)). Those are the next things to do in order to assess the validity of the current 
model and both of these goals is a natural continuation of the line of present studies. We will 
elaborate on observational properties and the implications of the time-dependent behavior 
of our solutions in the following paper of this series. 
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9. Conclusions 

Our studies of AGN infrared- driven winds have proceeded incrementally. In Paper I 
we considered a toy model which roughly approximated an infrared driven obscuring wind 
in AGN. An important prediction which resulted was that if the temperature of the dusty 
plasma in the rotating torus exceeds some characteristic value, the IR pressure on dust 
will inevitably overwhelm gravity, creating an outflow. However, this effective temperature, 
which approximately reads as T eff = (^) 1/4 ^ 312 (naMr)V4 _g 87 (nzMz)i/4 K> was derived 
analytically from a simplified model. 

In Paper I we combined calculations of stationary ID motion in the z-direction with 
2D calculations of the radiation field. The latter was calculated adopting 2D flux-limited 
diffusion approximation. In Paper II we relaxed the assumption of ID stationary motion 
and adopted a 2.5D time-dependent picture. One important "toy-model" simplification still 
remained though: there was no X-ray radiation explicitly taken into account. Instead, we 
calculated the effective temperature of the gas at the boundary of the computational domain. 
In the present work we relaxed this assumption and included X-rays into the global picture 
explicitly. 

Despite these advances, our work still does not consider the following potentially im- 
portant physical processes: 

The assumption of a razor thin disk is idealized. An accretion disk at parsec scales is 



likely self gravitating, possibly affected by gravitothermal instability (i.e., Gammie 2001). 



This may provide a local source of energy and maintain a Toomre parameter fix ~ 1 Toomre 



(1964). In a marginally gravitationally unstable thin disk, the non-axisymmetric perturba- 
tions lead to angular momentum transport which can be described by an effective viscosity 

Self-consistent models 



i.e., 



Shlosman et al. 


1989 


1990; 


Rafikov 


2009 



of geometrically thick disks in AGN at parsec scales call for the effects of self-gravity and 
radiation physics in 3D. This difficult endeavor can be a subject of future research. Another 
complication is that the problem of parsec scale winds may be intrinsically connected with 
the problem of accretion at parsec scale. Relaxing the assumption of an equatorial thin disk 
and treating accretion adopting an effective a prescription is a possible next step in the 
development of the current model. 

The assumption of a razor thin disk as a source of matter underpins the current study. 
Certainly such a distribution of matter is in many respects most unfavorable for the formation 
of radiation driven flow, owing to the geometry of X-ray illumination. Thus, we expect that 
if the wind does form for certain T in the thin disk case, then we may expect that in less 
idealized circumstances it will also form. To relax the assumption of a razor thin disk one 
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would need to abandon the assumption of equatorial symmetry and to include accretion disk 
into self-consistent consideration. 

Starting from the Paper I we see that much of the massive wind doesn't leave the 
potential well of the BH. This conclusion is confirmed through Paper II, and in the present 
work as well. The precise fate of this gas is beyond the scope of these works, but the results 
suggest the following: if the accreted gas is rejected by the BH at least some of it will return 
for additional attempts. 

We conclude with a mix of results and ideas driven by our current studies: 

1. The distribution of radiation force inside the dusty component of the flow depends 
on the shape of the photosphere where X-rays are converted into IR and can only be 
calculated in multi-dimensional simulations. This work includes both X-rays and IR 
into 2.5D time-dependent, hydrodynamic simulations . 

2. We do not observe a rotationally supported quasi-static torus in our simulations, al- 
though our simulated winds are likely to resemble some of the observed properties of 
such tori. 

3. The flow in the current simulations is more complex and time-dependent than that of 
our previous studies where we did not include X-rays explicitly. 

4. If geometrically thick dusty obscuration develops then a hot photo-ionized flow with 
velocities of 100 — 700 km s -1 accompanies it. X-ray warm absorbers are the evapo- 
rative flow originating both from the disk and also evaporated from such cold dusty 
component. 

5. We speculate that large scale motions with v z = ±(a few) x 100 km s" 1 originating 
from a spatially varying wind such we have calculated may be seen in maser emission 
in nearby bright type II AGNs. 

This research was supported by an appointment at the NASA Goddard Space Flight 
Center, administered by CRESST/UMD through a contract with NASA, and by grants from 
the NASA Astrophysics Theory Program 10-ATP10-0171. G.B-K. acknowledges the support 
of from the Russian Foundation for Basic Research (RFBR grant 11-02-00602). 
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10. Appendix A 

Frequency-integrated moments E, F, which appear in the above set of equations ([l])-((4]) 
are obtained by calculating angular moments from the frequency-integrated specific intensity, 
J(r, f2, u, t): 



E(r,t) = du (f) dQI(r,Q,u,t), (15) 

oo 



F(r,t) = J du j> dQhI(r,Q,u,t). (16) 
The frequency-independent radiation pressure tensor P is found from 

P(r, t) = ^J du j> dtt nn I(r, tt, v, t). (17) 
11. Appendix B 

Dust directly reprocess X-rays to IR and we approximately take this into account. The 
amount of energy absorbed by dust in a volume d 3 x: 

1 dE v f 

— = an d / a x I ux dQdu ~ an d a x cu x , (18) 



d 3 x dt 



where w x is the energy density of X-rays, cr x = x x /^d, and the number density of dust grains 
reads: rid = fdp/m>d, where m& is a mass of a grain, and fd is a dust-to-gas mass ratio. We 
make further simplifying assumption that the dust grain is being instantaneously heated to 
the effective temperature T d x . 

The contribution from the dust to the IR energy density (or, equivalently, to the tem- 
perature, T r is calculated from the following equation: 

= n d ca d (T djX - T r ) . (19) 
Assuming, that Td )X is constant over the time interval [t, t + dt], the update for radiation is 



T drX (t + dt) = T djX (l - e - dt / td ) + -e- d */' d T djX (t). Notice, that E = aT T 4 . The equation Q 
is used to in operator-split fashion to update E prior to all other terms in equation Q are 
updated. 



